char *tpmc_version = "!HVER!tpmc:   3.4.1 [CUED 12/22/11]";
char *tpmc_vc_id = "$Id: tmpc.c,v 1.1.1.1 2006/12/22 09:55:01 Exp $";

#include "HShell.h"     /* HMM ToolKit Modules */
#include "HMem.h"
#include "HMath.h"
#include "HSigP.h"
#include "HAudio.h"
#include "HWave.h"
#include "HVQ.h"
#include "HParm.h"
#include "HLabel.h"
#include "HModel.h"
#include "HTrain.h"
#include "HUtil.h"
#include "HAdapt.h"
#include "HFB.h"
#include "HDict.h"
#include "HNet.h"
#include "HLM.h"
#include "HLat.h"
#include "HArc.h"
#include "HFBLat.h"
#include "HExactMPE.h"
#include "exmath.h"
#include "xtpmc.h"
/* Parameters Global storage  */
static int numChans = 26;
static int numcepLifter = 22;
static int numceps = 12;
static MemHeap hmmNoiseStack;
static HMMSet hsetNoise;
static char *newDir = "";
/* Prototype function*/

void baoUpdateModels(DVector *mean, DVector *var, int numFeatures);

/* The following is HTK standard function*/
void Summary(void)
{
   printf("\nTPMC Command Summary\n\n");

   Exit(0);
}

void ReportUsage(void)
{
	printf("\nUSAGE: TPMC [options] Cleanpath Noisepath alpha ... \n\n");
	printf(" Option                                       Default\n\n");
	printf(" -mmfnoise  s path of noise model.\n");
	printf(" -hmmNoiseList  s path to list phones of noise.\n");
	printf(" -alpha  i the amount of noise used to corruped speechs.\n");
	printf(" -option  s trajectory PMC - simple average.\n");
	printf(" -numPoints  i if option ends with '_r', and numPoints > 0: generate sample points so 	that later HMM retrain could be done.\n");
	printf(" -noiseAddOpt  s this option is inactive so just set to \"normal\".\n");
	printf(" -covarOpt  s set to 'diag'.\n");
	printf(" -numObsVectors  i num of observation vectors (set to 2) to be.\n");
	printf(" -wOpt  s set 'l1-wNormal' reported in the paper (here, windowL=1).\n");
	printf(" -isDebug  i 0 (no output), 1 (some summary), 2 (full debug).\n");
	PrintStdOpts("Q");
	printf("\n\n");
}
/*LoadMMF: get mean and var in hset*/
void loadMMF(HMMSet *hset, Boolean ldBinary,DVector *stateNoiseMeanData,DVector *stateNoiseVarData, char **arNoiseInfo)
{

	HLink hmm;
	HMMScanState hss;
	int size;
	MixPDF* mp;

	int numComponents = 0; int numStates = 0; 	int numHMMs = 0;

	int numFeatures = hset->vecSize;

	int numMix = hset->numMix;

	DVector vTemp = NULL;
	vTemp = CreateDVector(&gstack,numFeatures);
	int ii = 1;
	NewHMMScan(hset, &hss);
	do {
		numHMMs++;
		hmm = hss.hmm;

		int stateId = 1;
		while (GoNextState(&hss,TRUE)) {
			stateId++;

			while (GoNextStream(&hss,TRUE)) {
				size = hset->vecSize;

				if (hss.isCont){

					int numMixtures = 0;
					numStates++;
					while (GoNextMix(&hss,TRUE)) {
						numMixtures++;
						numComponents++;

						TouchV(hss.mp->mean);

						for(ii=1; ii<=numFeatures; ii++)
							vTemp[ii] = hss.mp->mean[ii];

						stateNoiseMeanData[numComponents] = CreateDVector(&gstack,numFeatures);
						CopyDVector(vTemp,stateNoiseMeanData[numComponents]);
						TouchV(hss.mp->cov.var);

						for(ii=1; ii<=numFeatures; ii++){

							vTemp[ii] = hss.mp->cov.var[ii];
						}

						stateNoiseVarData[numComponents] = CreateDVector(&gstack,numFeatures);
						CopyDVector(vTemp,stateNoiseVarData[numComponents]);

					} /* while (GoNextMix(&hss,TRUE)) */

					char* phone = malloc(15*sizeof(char));
					char *strMixture = malloc(5*sizeof(char));;
					sprintf(phone, "h-%s-s%d", hss.mac->id->name, stateId);
					sprintf(strMixture, "%d", numMixtures);
					arNoiseInfo[3*numStates] = phone;
					arNoiseInfo[3*numStates+1] = strMixture;
					arNoiseInfo[3*numStates+2] = strMixture;

				} /* if (hss.isCont) */
			} /* while (GoNextStream(&hss,TRUE)) */
		} /*  while (GoNextState(&hss,TRUE)) */


	} while (GoNextHMM(&hss));
	EndHMMScan(&hss);
	/* load mean and var from tied*/

}
static HMMSet hset;                  /* Set of HMMs to be re-estimated */
static MemHeap hmmStack;           /* HMM defs and related structures */

int main(int argc, char *argv[])
{
	char *cleaninfo,*noisePrefix,*hmmNoiseList;
	int numPoints, numObsVectors, isDebug;
	double alpha;
	char *option, *noiseAddOpt, *covarOpt, *wOpt, *noisyPrefix;
	DMatrix dctDMatrix,cepDMatrix,invDctDMatrix, invCepDMatrix, Sq, W, VS;
	char*s;
	setbuf(stdout, NULL); /*unbuffered output.*/
	InitShell(argc,argv,tpmc_version,tpmc_vc_id);
	InitMem();
	InitMath();
	InitSigP();
	InitAudio();
	InitWave();
	InitVQ();
	InitModel();
	InitParm();
	InitLabel();
	InitTrain();
	InitUtil();
	InitFBLat();
	InitArc();
	InitDict();
	InitLat();
	InitNet();
	if (!InfoPrinted() && NumArgs() == 0)
	      ReportUsage();
	if (NumArgs() == 0) Exit(0);

	CreateHeap(&hmmStack,"HmmStore", MSTAK, 1, 1.0, 50000, 500000);
	CreateHMMSet(&hset,&hmmStack,TRUE);
	CreateHeap(&hmmNoiseStack,"HmmStore", MSTAK, 1, 1.0, 50000, 500000);
	CreateHMMSet(&hsetNoise,&hmmNoiseStack,TRUE);

	while (NextArg() == SWITCHARG) {
		s = GetSwtArg();
			switch(s[0]){
				case 'a':
					if(!strcmp(s, "alpha")){
						alpha = GetChkedFlt(0,10,s);

					}else HError(1, "Unknown option %s",s);

					break;

				case 'o':

					if(!strcmp(s, "option")) {
						option=GetStrArg();
					}
					break;

				case 'c':
								if(!strcmp(s,"covarOpt")){
									covarOpt = GetStrArg();
								}

								break;

				case 'h':
								if(!strcmp(s,"hmmNoiseList")){
									hmmNoiseList = GetStrArg();

								}
								break;
				case 'H':
									if (NextArg() != STRINGARG)
										HError(2319,"tpmc: HMM macro file name expected");
									if(!strcmp(s,"H")){

										char *x = GetStrArg();
										AddMMF(&hset,x);

									}
									break;
				case 'i':
									if (!strcmp(s, "isDebug")) {
										isDebug =  GetChkedInt(0,1,s);
									}
									break;
				case 'n':
									if (!strcmp(s, "noiseAddOpt")) {
										noiseAddOpt = GetStrArg();
									}else if(!strcmp(s, "numPoints")){
										numPoints = GetChkedInt(0,100,s);
									}else if(!strcmp(s, "numObsVectors")){
										numObsVectors = GetChkedInt(0,100,s);
									} else HError(1, "Unknown option %s", s);

									break;
				case 'm':
									if(!strcmp(s,"mmfnoise")){
										char *x = GetStrArg();
										AddMMF(&hsetNoise,x);
									}
									break;
				case 's':
									if(!strcmp(s,"savePath")){
										newDir = GetStrArg();

									}
									break;



				case 'w':
								if(!strcmp(s,"wOpt")){

									wOpt = GetStrArg();
								}

									break;

				default:
						HError(2319,"tpmc: Unknown switch %s",s);
						break;
			}
		} /*matches while(NextArg() == SWITCHARG)*/

	char *hmmList = GetStrArg();

	MakeHMMSet( &hset, hmmList );
	LoadHMMSet( &hset,NULL,NULL);
	SetParmHMMSet(&hset);

	MakeHMMSet( &hsetNoise, hmmNoiseList );
	LoadHMMSet( &hsetNoise,NULL,NULL);
	SetParmHMMSet(&hsetNoise);

	int numStates = hset.numStates;
	int numMix = hset.numMix;
	int vSize = hset.vecSize;
	int numStatesNoise = hset.numStates;
	int numMixNoise = hsetNoise.numMix;


	int M = numceps + 1; /* num statics features*/
	int numPhones = numStates;
	int totalNumFeatures = vSize;
	int D = totalNumFeatures / M; /* 1,...,(D-1)-order dynamic features*/
	int numDctDVectors;
	int windowL;

	DVector *stateCleanMeanData = malloc((numMix+1)*sizeof(DVector));
	DVector *stateCleanVarData = malloc((numMix+1)*sizeof(DVector));
	char **arCleanInfo = malloc((3*numStates+3)*sizeof(char*));

	DVector *stateNoiseMeanData = malloc((numMixNoise+1)*sizeof(DVector));
	DVector *stateNoiseVarData = malloc((numMixNoise+1)*sizeof(DVector));
	char **arNoiseInfo = malloc((3*numStatesNoise+3)*sizeof(char*));

	/*
	%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%% Preprocess options %%%
	%%%%%%%%%%%%%%%%%%%%%%%%%%
	*/

	char *origOption = option;

	processOptions(&option, &wOpt,&numPoints, numObsVectors,&windowL,&numDctDVectors);

/*	%%%%%%%%%%%%%%%%%%%%
	%% Initialization %%
	%%%%%%%%%%%%%%%%%%%%
*/

	DMatrix *Matrices,selMat;

	Matrices = initializeMatrices(numObsVectors, numDctDVectors, numChans, numceps, numcepLifter, M, D, windowL,option, wOpt);

	dctDMatrix = Matrices[1];
	cepDMatrix = Matrices[2];
	invDctDMatrix = Matrices[3];
	invCepDMatrix = Matrices[4];
	Sq = Matrices[5];
	W = Matrices[6];
	VS = Matrices[7];
	selMat = ones(1, numObsVectors);/* % duplicate numObsVectors time for other options except tpmc. For tpmc, we will reset it later*/

	puts("Pre-process Noise statistics... ");
	loadMMF(&hsetNoise, 0,stateNoiseMeanData,stateNoiseVarData,arNoiseInfo);


	DMatrix noiseMeanLinear, noiseVarLinear, noiseMeanLog, noiseVarLog, noiseMeanMfcc, noiseVarMfcc, trajMeanVect, trajVarMat;
	DMatrix *preProDMatrix;

	noiseMeanMfcc = CreateDMatrix(&gstack,1,totalNumFeatures);
	noiseVarMfcc = CreateDMatrix(&gstack,1,totalNumFeatures);
	CopyDVector(stateNoiseMeanData[1],noiseMeanMfcc[1]);
	CopyDVector(stateNoiseVarData[1],noiseVarMfcc[1]);

	preProDMatrix = preprocess(noiseMeanMfcc[1], noiseVarMfcc[1], M, Sq, numDctDVectors, D, W, windowL, invDctDMatrix, invCepDMatrix, covarOpt, option, selMat, isDebug,0);
	noiseMeanLinear = preProDMatrix[1];
	noiseVarLinear = preProDMatrix[2];
	noiseMeanLog = preProDMatrix[3];
	noiseVarLog = preProDMatrix[4];
	noiseMeanMfcc = preProDMatrix[5];
	noiseVarMfcc = preProDMatrix[6];
	trajMeanVect = preProDMatrix[7];
	trajVarMat = preProDMatrix[8];

	/*%%% Main process %%%*/
	int startId,debugStartId,debugEndId,prevNumMixture,orignumObsVectors,phoneId,numDefaultMixtures,numMixture,endId;
	char *phone;
	startId = 1;
	debugStartId = 1;
	debugEndId = 1;
	prevNumMixture = -1;
	orignumObsVectors = numObsVectors;
	DMatrix noisyMeanData,noisyVarData;

	noisyMeanData = CreateDMatrix(&gstack,numMix,M*D);
	noisyVarData = CreateDMatrix(&gstack,numMix,M*D);

	loadMMF(&hset, 0,stateCleanMeanData,stateCleanVarData,arCleanInfo);

	int k;
	DMatrix *stateNoisyData;
	int i,j;
	DMatrix factDMatrix;/*=kron(ones(1, numObsVectors), eye(numFeatures,numFeatures));*/

	factDMatrix= CreateDMatrix(&gstack,M*D,M*D*numObsVectors);
	for(i=1;i<=M*D;i++)
		for(j=1;j<=M*D*numObsVectors;j++)
			if(i==j || M*D+i ==j )
				factDMatrix[i][j]=1;
	char *val=malloc(30*sizeof(char));

	clock_t start = clock();
	DVector *CleanMeanData = malloc(32*sizeof(DVector));
	DVector *CleanVarData = malloc(32*sizeof(DVector));

	/*ShowDVector("",stateCleanMeanData[101],39);*/
	/*ShowDVector("",stateCleanVarData[101],39);*/
	DVector memMarker;
	puts("Start TPMC converting...");
	for(phoneId=1;phoneId<=numPhones;phoneId++)
	{

		phone = arCleanInfo[3*phoneId];
		numDefaultMixtures = atoi(arCleanInfo[3*phoneId+1]);
		numMixture = atoi(arCleanInfo[3*phoneId+2]);
		endId = startId + numMixture - 1;
		for(i=0;i<numMixture;i++)
		{
			CleanMeanData[i+1]= stateCleanMeanData[startId+i];
			CleanVarData[i+1]= stateCleanVarData[startId+i];
		}

		if(phoneId%100==0)
		 {
			printf( "Done processing %d out of total %d states...\n",phoneId,numPhones);
			fflush(stdout);
		 }

		memMarker=CreateDVector(&gstack,1);
		stateNoisyData = modelPMC(phoneId, phone, noiseMeanLinear, noiseVarLinear, noiseMeanLog, noiseVarLog, CleanMeanData, CleanVarData, numMixture, alpha, M, numObsVectors, numDctDVectors, D, W, VS, Sq, windowL, numChans, numceps, dctDMatrix, cepDMatrix, invDctDMatrix, invCepDMatrix, covarOpt, option, selMat,factDMatrix, isDebug);


		for(k=1;k<=numMixture;k++)
		{
			CopyDVector(stateNoisyData[1][k],noisyMeanData[startId+k-1]);
			CopyDVector(stateNoisyData[2][k],noisyVarData[startId+k-1]);
		}

		FreeDVector(&gstack,memMarker);

		/*% update variables*/
		startId = endId + 1;

	}/*end phone*/


	printf("Done TPMC using %f min\n", ((double)clock()-start)/(60*CLOCKS_PER_SEC));
	fflush(stdout);
	FreeDMatrix(&gstack,factDMatrix);

	/*ShowDVector("",noisyMeanData[101],39);*/
	/*ShowDVector("",noisyVarData[101],39);*/
	puts("Save noisy stats to MMF format");
	baoUpdateModels(noisyMeanData,noisyVarData,totalNumFeatures);

	return (0);          /* keep compiler happy */
}
void baoUpdateModels(DVector *mean, DVector *var, int numFeatures)
{
   HMMScanState hss;
/*cpdf = pdf[s].spdf.cpdf = CreateCME(hset,M);*/
   int i;
   {
      int k=0;
      if(hset.hsKind==PLAINHS || hset.hsKind==SHAREDHS){
         NewHMMScan(&hset, &hss);
         while(GoNextMix(&hss, FALSE)){ /* Update mixes. */
            k++;
            for(i=1;i<=numFeatures;i++)
            {
            	hss.mp->mean[i] = mean[k][i];
            	hss.mp->cov.var[i] = var[k][i];
            }
            /*FixGConsts(hss.hmm);*/
            FixDiagGConst(hss.mp);
         }
         EndHMMScan(&hss);
      }
   }

   printf("Saving hmm's to dir %s\n",(newDir==NULL)?"Current":newDir);
   fflush(stdout);
   SaveHMMSet(&hset,newDir,NULL,NULL,0);
}








